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We consider the set of Diophantine equations that arise in the context of the partial 
differential equation called "barotropic vorticity equation" on periodic domains, when 
nonlinear wave interactions are studied to leading order in the amplitudes. The solutions 
to this set of Diophantine equations are of interest in atmosphere (Rossby waves) and 
Tokamak plasmas (drift waves), because they provide the values of the spectral wavevec- 
tors that interact resonantly via three-wave interactions. These wavenumbers come in 
"triads", i.e., groups of three wavevectors. 

We provide the full solution to the Diophantine equations in the physically sensible 
limit when the Rossby deformation radius is infinite. The method is completely new, 
and relies on mapping the unknown variables via rational transformations, first to ratio- 
nal points on elliptic curves and surfaces, and from there to rational points on quadratic 
forms of "Minkowski" type (such as the familiar space-time in special relativity). Classical 
methods invented centuries ago by Fermat, Euler, Lagrange, Minkowski, are used to clas- 
sify all solutions to our original Diophantine equations, thus providing a computational 
method to generate numerically all the resonant triads in the system. Computationally 
speaking, our method has a clear advantage over brute-force numerical search: on a 
10000 2 grid, the brute-force search would take 15 years using optimised C ++ codes on 
a cluster, whereas our method takes about 40 minutes using a laptop. 

Moreover, the method is extended to generate so-called quasi-resonant triads, which 
are defined by relaxing the resonant condition on the frequencies, allowing for a small 
mismatch. Quasi-resonant triads' distribution in wavevector space is robust with respect 
to physical perturbations, unlike resonant triads' distribution. Therefore, the extended 
method is really valuable in practical terms. We show that the set of quasi-resonant 
triads form an intricate network of connected triads, forming clusters whose structure 
depends on the value of the allowed mismatch. It is believed that understanding this 
network is absolutely relevant to understanding turbulence. We provide some quantitative 
comparison between the clusters' structure and the onset of fully nonlinear turbulent 
regime in the barotropic vorticity equation, and we provide perspectives for new research. 
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I. INTRODUCTION 

Understanding the behaviour of nonlinear wave resonances is of importance in numerical 
weather prediction, fusion reactors and in general in experiments in fibre optics and water waves. 
Modulational instability is the first example of relevance of resonant and quasi-resonant triads 
and quartets [lj. Effects such as wave turbulence with its cascades of energy and enstrophy rely 
on the interconnected web of resonances. However a thorough understanding of the full-scale 
mechanism of energy transfers is lacking at the moment and one of the reasons for this is the 
poor understanding of the web of resonant and quasi-resonant triads. This poor understanding 
is in two aspects: (i) Kinematical aspect, regarding the actual position in wavenumber space of 
the modes that interact in resonant and quasi-resonant triads, as well as the interconnections of 
these triads to form clusters, (ii) Dynamical aspect, regarding the processes of energy exchanges 
between modes via 3-wave interactions. This paper will concentrate on the kinematical aspect, in 
the context of the barotropic vorticity equation on periodic domains. As a motivation, we believe 
that the information obtained from the kinematical aspects can be used to derive a reduced 
approximated system of evolution equations for the Fourier amplitudes of the underlying PDEs 
that describe the nonlinear oscillations. The long-term goal is to produce an important degree 
of simplification from an infinite-dimensional space of dynamical variables to a finite-dimensional 
one, while keeping enough complexity so as to accommodate physically-measurable characteristics 
such as turbulence and energy/enstrophy cascades. This paradigm is in the spirit of the study of 
mesoscopic and discrete wave turbulence (areas under development). See for example [5HS] . 

Another question is whether energy transfers between wave modes depend or not on the 
nature and strength of the forcing. Such question has been studied in a myriad of papers using 
numerical simulations of full PDE models, and the intuition gained is quite empirical. However, 
the understanding of the network of connected triads might lead to a better modelling of this 
and other problems. In [6], the effect of forcing a single triad was studied and it was established 
that the triad's energy remains bounded, even if the forcing oscillates in phase with the natural 
frequency of the so-called unstable mode. In [7J, formation of Kolmogorov-like spectra of inverse 
cascade of weak turbulence was observed in a direct numerical simulation of a special formulation 
of water waves. This formulation, known as integro-differential Zakharov equation, was exploited 
numerically in the cited reference so that only a relatively small set of wave modes needed to be 
integrated forward in time. 

The importance of the quasi-resonant triads over the exactly resonant ones for the accurate 
and faithful description of nonlinear physical wave systems has been verified in several works, 
starting from the modulational instability [lj and nowadays it is an accepted fact that any proper 
description of wave interactions must take into account quasi-resonant triads from the start 

mini. 

In Jupiter's atmosphere, local wave oscillations have been observed with wavelengths of 20 
km [? ]. If such waves could be extended to cover a significant part of Jupiter's surface then 
we would have wavenumbers of approximately 3000. In numerical simulations, particularly in 
atmospheric models, it has been shown [TT| that high resolutions are needed if one wants to 
describe properly the generation of small-scale turbulence out of a large-scale initial condition 
(see the benchmark problem of Rossby-Haurwitz wave in the cited reference). Therefore, running 
2D simulations with up to 10000 grid points per spatial direction is meaningful, particularly if 
one is interested in studying inverse energy cascades. The problem with these high-resolution 
simulations, is that a rigorous analysis of the wave interactions (particularly in the limit of very 
small amplitudes, used in wave turbulence) escapes from the computational capabilities of the 
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current computer technology. For example, if we try to do a numerical search of the quasi- 
resonant modes with wavenumbers less than or equal to 10000 in size, a computer will need to 
do a search over 10000 4 possible triads, and evaluate the frequency mismatch. If r is the time 
it takes the computer to do one triad, then the total time for this search will be approximately 
r x 10 16 . On a medium-sized cluster, using optimised C ++ codes, the effective time for one such 
computation is about rs;5x 10~ 8 seconds, which means that the total computing time for all 
triads is about 15 years! 

The search for quasi-resonances becomes therefore a computational problem and this paper 
provides a solution, with a computational method that is based on a mapping from resonant 
triads to rational points on quadratic hyper-surfaces of "Minkowski" type (i.e., non-definite). The 
method for searching exact resonances is extended to a method for finding quasi-resonant triads. 
Here the advantage of the method over brute-force direct search is that the method produces 
quasi-resonant triads with preferentially small detuning levels, which is useful physically, since we 
are interested in the regime of small enough amplitudes. A Mathematica code that generates 
the triads using this method is available from the authors upon request. 

The paper is organised as follows. In Section [M] we review the barotropic vorticity equation 
and the Rossby triads on a periodic spatial domain (torus). In Section III we state the Diophan- 
tine equation that determines the triad resonance conditions, develop the mapping to elliptic 
curves and to Minkowski-type of quadratic forms, establish the general solution to the resonance 
conditions, and explain in detail the computational method of solution, with a step-by-step for- 
mulation. In Section ]V we present the numerical application of the method, showing how to 
construct the triads in the box of size 5000 in wavenumber space, and studying the structure, 
distribution and connectivity properties of the resonant clusters that form. In Section |V]we in- 
troduce quasi-resonant triads, which are far more relevant physically than resonant triads. We 
explain the method to construct these quasi-resonant triads, which is based on the previous 
method of construction of resonant triads. Finally, in Section VI we discuss further extensions of 
the method to the case when the aspect ratio of the spatial variables is not equal to one. The 
Appendix contains all necessary mathematical results to understand the method of construction 
of resonant triads presented in Section pTTj 



II. GOVERNING EQUATIONS 
A. Rossby Triads 

We consider the dynamics of a rotating shallow layer of incompressible fluid, in the so-called 
/3-plane approximation. The governing equation reduces to the following partial differential 
equation, known as the barotropic vorticity equation: 

| ( VV-W§^-^)+4^0, (!) 
at v ' \ ax oy ay ox J ox 

where ip = -ip(x,y,t) 6 K is the streamf unction, (3 is a constant that determines the speed of 
rotation of the system, and F = 1/R 2 where R is the Rossby deformation radius. For simplicity, 
we consider from here on the physically sensible limit of infinite Rossby deformation radius, 

F = 0. 

Equation ([I]) is also known in the literature as the Charney-Hasegawa-Mima equation, or 
CHM. The "C" comes from the atmospherical context just described. The "HM" comes from 
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an independent derivation of the same mathematical equation in the context of plasma physics 
[12]. Due to this multidisciplinary aspect, the CHM equation occupies a special place in physics 
and mathematics. On the one hand, it suffices to say that there is no general solution of the 
equation. Due to its nonlinear term, its solutions display an amazing variety of complex behaviour 
and turbulence, such as inverse energy cascades (from small to large scales), direct enstrophy 
cascades (from large to small scales) and zonostrophy cascades (with non-isotropic transfers). On 
the other hand, the model is simple enough to accommodate convenient mathematical features 
such as a non-canonincal Hamiltonian structure [13], and the equation conserves total energy and 
enstrophy. More importantly for our study of resonances, the equation admits an infinite number 
of so-called "travelling wave" solutions. These are simply solutions of the linearised version of ([l]) 
in the form of a travelling wave (called Rossby wave), parameterised by wavevectors (k,l) £ M : 

*(*,,)(z,y,i) = e i(lI+ ' s -" (i ' 0,> , (2) 

«(*.0— j^p- (3) 

It is easy to check that ^(k,i), as well as its real part, satisfies also the full equation ([I]), i.e., the 
nonlinear term in that equation is identically zero. However, an arbitrary linear cornEination of 
travelling wave solutions with wavevectors that are not collinear, is still a solution of the linearised 
version of ([I]) but is not anymore a solution of the full equation. 

The last observation leads to the study of the nonlinear term and the mixing of modes, and is 
the starting point in the construction of so-called resonant triad solutions, which are approximate 
solutions, valid in the asymptotic limit when the oscillation amplitudes are small. Without going 
into the details of the multiple-scales method [14], these solutions are linear combinations of 
three travelling waves, of the form: 

&(x,y,t) = &(A 1 (t)V (kl _ M (x,y,t) +A 2 (t)y {k2th) (x,y,t) + A 3 (t)* {ka , la) (x,y,t)) , (4) 

where 3? denotes real part, A, are some complex functions of time that are "slow" compared to 
the waves, and the set of wave vectors satisfy the following system of equations: 

h + h = k 3 (5) 

h + k = h (6) 
wi + uj 2 = uj 3 , (7) 

where Uj = u(kj, lj) , j — 1, 2, 3 . Any set of three wavevectors satisfying equations (j5])-(|7]) is 
called a resonant triad. We remark that there is a consistency condition that leads to a nonlinear 
system of evolution equations for the three slow functions A,. We refer the reader to the vast 
literature on the subject, such as (THE] and references therein. In a more general setting, these 
slow functions also depend on the space coordinates [16], again 'more slowly' than the travelling 
waves. 

It is possible to extend the resonant triad solutions to a more general approximate solution, 
consisting of a combination of groups of resonant triads: 

N 

®(x,y,t) = $tJ2Mt)V(k I ,i I )(x,y,t) , (8) 

7=1 



in such a way that any wavevector (ki, h) appearing in the combination above, belongs to at least 
one resonant triad. In this setting, several resonant triads may be connected, forming so-called 



5 



clusters. The main question of this paper is the classification of all possible values of wavevectors 
(kj,lj) that belong to resonant and near-resonant triads when the spatial domain is periodic. 
However, it is worth mentioning that the N functions Aj(t) satisfy a coupled nonlinear system 
of evolution equations, that is derived using a set of consistency conditions, in a straightforward 
manner [171 US]. 



B. Periodic Domains and Numerical Simulations 



Periodic spatial domains are the common place for numerical simulations. One of the most 
reliable existing methods for periodic domains is the pseudo-spectral method, that takes ad- 
vantage of the fast Fourier transform to compute partial derivatives and nonlinear terms. For 
quadratic nonlinearities, as in the CHM equation, conservation of energy as well as enstrophy 
is guaranteed. The pseudo-spectral method is widely used in numerical experiments due to its 
rapid convergence in terms of accuracy, and is preferred over finite-difference and finite-volume 
methods for fundamental studies, such as turbulence and weak wave turbulence. 

In the schematic version of the pseudo-spectral method, solutions to equation ((!]) are sought 
on the periodic domain (x,y) E [0,27r) x [0, 2n), so solutions ip(x,y,t) satisfy ip(x,y,t) = 
ip(x + 2n,y,t) = ip(x,y + 2ir,t). The numerical method begins by approximating the solution 
ip(x,y,t) as a finite sum of Fourier modes: 

N X Ny 

i[>(x,y,t) =$Y1 Yl A (*.')(*)*(fc,!)O c >y)*)) ( 9 ) 

k=0 l=-N y 

where a "mode" ^^,i)(x,y,t) is defined in equation ([2]) and is identified with its associated 
wavevector (k,l). The unknown coefficients A/f.,i){t) nee d to be advanced in time numerically, 
using high-order techniques such as 4* ft -order Runge-Kutta. The natural numbers N x ,N y denote 
the spatial resolution of the numerical approximation, so that oscillations of wavenumber greater 
than these numbers cannot be resolved by the numerical scheme. 



III. DISCRETE RESONANT TRIADS 



In the periodic setting discussed in Section MB the relevant wavevectors (k,l) belong to 
a discrete lattice of integer numbers, so the resonant conditions (p])-(Kl) are now a system of 
Diophantine equations. For Rossby waves, the dispersion relation (|3]) applies, and after some 
algebra the resonant conditions can be reduced to the following equation: 



HK + ify - 2k 3 (ki + if)(k 3 h + 1 3 h) + 2h(k 3 h + 1 3 h)(ki + ii) - h(ki + liy = o , (io) 

which can be interpreted as an equation for the pair (ki,li), if the pair (k 3 ,l 3 ) is given. Here, 
the pair (A^,^) ' s obtained a posteriori by writing (k2,h) — (^3,^3) — (ki,h). 



Discarding Zonal Modes. A mode (k,l) with k 
k 3 = 0, equation (10) reduces to (21 1 — h)ki = 



= represents a so-called zonal mode. When 
0, which is easy to solve. This case is not 
considered here because resonant interactions with zonal modes are trivial: physically these 
interactions are suppressed due to the vanishing of interaction coefficients, see for example [TH1 
equation 13.9]. We remark that zonal flows are interesting, though, if quasi-resonant interactions 
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are allowed (see Section [V]). 

So, from here on we will assume ki,k 2 ,k 3 ^ 0. Without loss of generality one can assume 
also that < k\ < k 2 < k 3 , but this latter assumption is not essential until the counting of 



physically different solutions of equation (10) is performed 



Irreducible Triads. Notice that equations (10) are invariant under overall re-scaling of the 



wavenumbers. This is useful from the computational point of view, because it means that we 
need to search for so-called irreducible triads, which are defined as solutions of the resonant 



condition (10) such that the six wavenumbers kx, l\, k 2 , h, k 3 , ^3 do not have a common factor. 
Correspondingly, a reducible triad is defined by any solution that is an integer multiple of an 
irreducible triad. 

Change of Basis. Let us make a convenient change of basis to the orthogonal basis spanned 

by (k 3 ,l 3 ) and (-l 3 ,k 3 ): 

(k 1 ,h)=a(k 3 ,l 3 ) + b(-l 3 ,k 3 ), (11) 
where (a, b) E Q 2 . Explicitly, the inverse transformation is 

hh + hh k 3 h-l 3 h 
a = k 2 + p > b = k 2 + p ■ ( 12 ) 

3 3 3 3 



Notice that kf + lj = (a 2 + b 2 )(kl + ll) and k\ = ak 3 -bl 3 . Therefore, equation ( |l0[ ) transforms 
to: 

{a 2 + b 2 ) 2 - 2a(a 2 + b 2 ) + (2a - 1) (a - ^ bj = . (13) 

This equation is then equivalent to the original resonant condition, provided zonal modes are not 
involved. 

A. Mapping the resonance condition fll3| ) to classical problems: pure-cube solutions and 

elliptic curves 

Two cases need to be discussed separately: 



Case a = 0. Equation (13) simplifies to 



6 4 + ii b = . (14) 

k: 



3 

We discard the solution 6 = because it represents a triad that is formed by collinear modes, 
and in this case the interaction coefficients are identically zero (i.e., the mixed contribution 
stemming from the nonlinear terms in equation ((!]) vanish identically). So we take, without loss 
of generality, 6^0. We get then b 3 = — j^. But b is a rational number, so we deduce that the 
case a = gives a nontrivial triad if and only if the ratio ^ is a pure-cube rational. In such 
instance we obtain a triad of the form (ki, h), (k 2 , l 2 ), (k 3 , l 3 ), where 

/ \ 1/3 



(A*, h) = -[-±] (-l 3 , k 3 ), (k 2 , l 2 ) = (k 3 , l 3 ) - (k u h), (15) 
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where an overall scaling factor might be needed in order that all wavenumbers be integer. 

/ \l/3 I I 

Example: Z 3 = l,k 3 = 8. Then b = — ( |M = — | and so, applying equation (15) directly 
we get a preliminary triad: 

(Kh)' = (\, -4) , (k 2 , h)' = , 5] , (* 8 , is)' = (8, 1) , 

so we multiply by the factor 2 in order to get an integer irreducible triad: 

(k 1} l x ) = (1, -8) , (k 2 , l 2 ) = (15, 10) , (k 3 , l 3 ) = (16, 2) . 

This "pure-cube" construction can be done for arbitrarily large values of wavenumbers and 
it leads to an infinite set of different irreducible resonant triads. But most of the solutions of 
the resonant conditions (13) are found via a completely different method, to be discussed in the 
remaining of this Section. 

Case a 7^ 0. Provided the ratio l 3 /k 3 is not a pure-cube rational, we define £ = - and thus 
equation (13) becomes 

S(1 + _ 2 Q 2 (1 + ^ + (2 a _ 1} A _ I ^ = (16) 

Notice that this equation has degree 3 in a, one less than the original equation. We follow a 
standard procedure to isolate the coefficient of the cubic term. Defining r = a(l + £ 2 ) we get 

,3 o 2 1 / „ ( 1 1 f2\\ / -1 ^3 



r 3_2r 2 + (2r-(l+£ 2 )) ( 1 - 1 = 0. (17) 

Defining now D = — and the quantities X = rD, Y = £D , we obtain the equation 

X 3 -2DX 2 + 2DX - D 2 = Y 2 , (18) 
which can be interpreted as an Elliptic Curve if D is fixed and we let X, Y be variable. 

Summary of case a 7^ 0. To summarize, provided the ratio h/k 3 is not a pure-cube rational, 
there is a rational mapping from an integer triad (ki,l%), (k 2 ,l 2 ), (^3,^3) satisfying all resonance 
conditions (|5])-([7]), to the variables of the elliptic curve (18). The mapping is bijective up to 
overall re-scaling of the triad wavenumbers, and given explicitly by: 

fc 3 k 2 + 1 2 v _h k 3 h - k x l 3 n _h hh + hh 



and with inverse 

fci X h X f D \ l 3 D - 1 



A; 3 £> 2 + y 2 ' fc 3 y V -D 2 + y 2 / ' fc 3 y 

and (k 2 ,l 2 ) = (k 3 - h,l 3 - h). 



(20) 
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B. Classification of solutions of triad equations in terms of Fermat's theorem of sums of 

squares 

We consider the case a ^ described in the previous Section. We can rewrite the elliptic 

(21) 



curve (18) as 



Y 2 + (D + X 2 - X) = X 2 (X 2 - X + l) , 
and we can divide by X because of our assumption of no zonal modes. We obtain 

2 / ^ \ 2 

X-l\ =X 2 -X + 1. 



D 

X 



(22) 



The LHS of equation (22) is a sum of squares of rationals. The RHS is a quadratic form that is 
best written in diagonal form by defining 

X = , m,neZ, (23) 



m 



n 



so we get the equation: 

Y(m 



n 



m + n 



+ 



Dim 



n) 



m + n 



2 m 



3 m 2 + n 2 



(24) 



Equation (24) is to be solved for m,n,£ Z and Y, D G Q. The cases Y = 0, m = or 
m = ±n are excluded from the solutions: they give rise to zonal modes. Notice that since 



m,n are integers, the two members of equation (24) are equal to an integer. The problem 
of finding all possible integers that can be written as a sum of squares of two rationals is well 
known and dates back to Fermat [20, Chapter V] (it is called Fermat's Xmas theorem). Also, 

was considered by Fermat and 



n 



the problem of finding all possible integers in the form 3m 
solved by Lagrange and Euler, see [21, Chapter 1] for further details. The current equation is 
a combination of these two problems and can be dealt with in a straightforward manner. As a 
result, all possible representations for the numbers m,n and the numbers Y,D can be obtained 
explicitly. Consequently, all possible exact resonant triads can be obtained explicitly by using the 



mapping (20) from X,Y,D to the triad's wavenumbers. 



Method for finding the general solution of equation (24). We believe that it is not 



illuminating to provide here a detailed exposition of the solution method. We have relegated to 
the Appendix all necessary Theorems and Corollaries. Here we provide only the main idea of the 
solution method, but we still give the detailed explicit construction of the solution. The main 



idea is that integers of the form (24), LHS, must be products of prime numbers of the form 
AK + 1 with K some integer, times a square. On the other hand, integers of the form (24), 
RHS, must be products of prime numbers of the form 3K' + 1 with K' some integer, times 



another square. So equating the two members of equation (24) we deduce that these must be 
equal to an integer that is product of primes of the form 12if" + 1 with K" some integer, times 
a square. 



It follows from Corollary 4 in the Appendix, that all possible solutions of equation (24) for 
m,n G Z coprime and Y, D G Q, can be parameterised explicitly using the following algorithm: 



(i) Construct the following expansion in prime powers: 



M 



N 



1 "0 



X 



Up7 x 



i=i 



' M' 

ft 



X 



'M" 

n 



7 j 



(25) 
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where: 

a can take the values or 1; 

{pj}JL l is the set of primes of the form 12k + 1 with k integer; 

{qj}jLi is the set of primes of the form 3 k + 1 with k integer, excluding the primes pf 

{ r j} C jLi is the set of primes of the form Ak + 1 with k integer, excluding the primes pf 

The non-negative integers M,M',M" denote how many primes are there in the expansion 
of N: the explicit expansion is given via three strings of non-negative integers {dj}jLi, 
{bj}jL 1 and {cj}jL v These are defined so that if M = then the string {aj}jL 1 is empty, 
and if M > then au > 0. Similar relations hold for M replaced by M' (resp. M") and 
cij replaced by bj (resp. Cj). 

Once a and the strings {aj}jL lt {bj}jL 1 and {cj}jCi a re known, the corresponding is 
uniquely constructed. 

(ii) Construct all possible solutions m, n e Z of the equation 



N = 3m 2 + n 2 



(26) 



using the Brahmagupta identity (Al) on the individual expansions of the solutions of 

3xl 2 + F 



Pj = 3 m 2 + n 2 - and 4 



The choice of sign for m leads to two sets of physically sensible solutions, and one can 
take n > without loss of generality. With this taken into account, it can be shown that 
the number of different solutions for m, n (with m^O, ±n) is equal to 



2 x (2a + 1) x - x 



' M M' 

HH(a j + l)(2b k + l) 
j=i k=i 



where e = 1 if all {dj}j =l are even, and e = otherwise, 
(iii) Construct all possible solutions S,Q G Z of the equation 

N = S 2 + Q 2 , 



(27) 



using the Brahmagupta identity (Al) on the individual expansions of the solutions of 

p 3 =S 2 + (.)]. 

We can take without loss of generality the convention < S < Q. With this taken into 
account, it can be shown that the number of different solutions for S, Q is equal to 



' M M" 

nn^ +i ) ( 2c '+ i )- e 
j=i i=i 



where e was defined in point (ii). Notice that the factor 4 a ° does not play an important 
role here. 
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(iv) Based on the above solutions for S, Q integers, construct all possible solutions s,q G Q\Z 
of the equation 

N = s 2 + q 2 , (28) 



using the Brahmagupta identity (Al) on the solutions obtained in point (iii) along with 



the Pythagorean rationals, which are defined by the remarkable identity: 

Vu,veZ. (29) 



. i \ - / 2 2 \ 2 

lUV \ I U — V x 



V? + V 2 J \U 2 + V 2 

The relevant solutions for s, q are of the form: 

2uvQ + (u 2 -v 2 )S -2uvS+ (u 2 - v 2 ) Q 

S = 2~7~2 ' 1 = 2~7~2 ' 30 ) 

where u, v are coprime integers satisfying 1 < |u| < u. This ordering prevents repetition 
of solutions, which at this stage we need to avoid, although they play a role in point (v). 

For computational purposes we parameterise the coprime integers u, v using finite sets of 
adjustable size. Letting U max G Z, U max > 2, define the "Pythagorean Fan" F(U max ) by 
the equation 

F(U max ) = {(u,v) G Z 2 coprime s.t. 1 < \v\ < u < U max } . 



Then, for a given solution S, Q of equation (27), a set of non-integer rational solutions s, q 



of equation (28) is obtained which amounts to #F(C/ max ) new solutions, where #F(U max ) 



is the cardinality of F(U n 



(v) Gather all results above and construct all possible solutions for Y and D of equation (24). 
Here, an extra 8-fold symmetry is available that allows us to generate more solutions and 
correspondingly more triads. Notice that for each obtained solution s,q (integer as well as 



non-integer) of equation (28), and solution m, n of equation (26), the following statement 
is true: Provided sq ^ 0, we have the choice of writing eight different equations for Y, D, 
schematically written as: 

Y(m — n) 2 Dim — n) 2 

— ~ = ±s, — - + 2m = ±q, 31 

m + n m + n 

and 

Y(m — n) 2 D(m — n) 2 

— ~ = ±q, — - + 2m = ±s, 32 

m + n m + n 

where all sign combinations can be taken. Notice that the symmetry Y —Y amounts 
to taking the mirror image of a triad about the /c-axis in the (k, Z)-wavevector space (see 



Eq. ©). 

There is a special case which occurs when all a/s are even. There, amongst all possible 
solutions S, Q there is one with S = 0. Moreover, even when SQ ^ there is a solution 
s,q with sq = for some choice of u, v. In these cases we still have a symmetry, but it 
reduces to a 2-fold symmetry, with two different equations for Y,D: 

Y(m — n) 2 Dim — n) 2 , , 

— ^ - = ±s, — J - + 2m = 0. 33 

m + n m + n 
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Gathering all solutions together, we can write a formula for the total number of solutions for 



X,Y,D of equation (21), leading to the following number of irreducible triads as a function of 



M f M' \ z / M" \ ' 

the integer N = 4 ao x Y[p a / x ( flff* ) x ( EM* J 

3=1 \k=l J \l=l J 



T(N) = (2a + 1) x 

M M" 



m W 

nn^ +i ) (2& fc +i)-e 

U=l k=l 



nn^ +i )( 2c ' +i )+ e 

j=i i=i 



x (4#F(t/ max )+4-3e). 



(34) 



For example, taking N = 4 x 13 and U max = 2, so a = l,ai = 1,M = 1,M' = M" = 
0, #F(f/ max ) = 2,e = 0, we get T(A r ) = 144 triads. 

Remarks. 

• As defined earlier, an irreducible triad is one whose wavenumbers ki,h,k 2 , h,k3,h do 
not have a common factor. Knowing all irreducible triads within a given box in (k,l)- 
wavenumber space, say \k\, \l\ < L, allows one to compute all reducible triads within that 
box, by including integer multiples of irreducible triads of size smaller than L/2. 



M 

Nested character of the triads. Let the integers Ni,N 2 be of the form 4 a ° x YIp^ 

3=1 

k=i / \i=i 

obtained using Ni contain the set of triads obtained using A^. 



x 



Wrf J . Then if A - ! = N 2 x N 2 for some integer N 3 , then the set of triads 



Eliminating extra multiplicity. The set of all possible triads (using all possible values of N) 
obtained from this method, will appear "repeated six times" in the sense of the six-fold 
permutation symmetry of the triad equations: 



ki + k 2 = k 3 , k 2 + kx = k 3 , 
-kg) +k 2 = (-ka), k 2 + (-k 3 ) = 
-kg) +ki = (-k 2 ), k x + (-k 3 ) = 



-ki) 
-k 2 ) 



(35) 
(36) 
(37) 



where k, = (kj,lj). The symmetry reduces to four-fold when permutations of "pure-cube" 
triads occur, because a "pure-cube" triad cannot be mapped to an elliptic curve, but four 
of its permutations can be mapped. Correspondingly, when collecting all triads obtained 



using our method, a normalisation must be applied if one wants to use formula (34) to 
estimate the number of genuinely different triads. The exact number of triads for a given 
A^ is thus a bit different than what the formula predicts. Only when triads are collected 



using several values of N, the formula (34) will give asymptotically a good estimation if 
we divide by 6. As an example, the choice A^ = 4 x 13 gives 40 genuinely different triads, 
rather than 144/6 = 24. 
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IV. NUMERICAL RESULTS: GENERATING RESONANT TRIADS WITH 
WAVEN UMBERS (k,l) IN THE BOX \k\ < 5000, |Z| < 5000 

Let us consider the following problem: Find all resonant triads in the box \k\, \l\ < 5000. As 
explained in the introduction, a direct search using a brute-force algorithm would take more than 

15 years to answer that problem. With the method introduced in this paper, such problem is 
reduced to the problem of finding the appropriate set of primes pj, qj, rj that generate all triads 
within the box. 

The answer to this latter problem has not been presented in this paper, and we believe this 
is a matter for a subsequent paper. From the rational character of the mapping from the triads 
to the points X, Y, D on the elliptic curve, it follows that only small-enough primes Pj,qj,rj 
and small-enough powers aj,bj,Cj can generate small resonant triads. The open question is how 
small should the primes and the powers be. An extra complication is given by the Pythagorean 
Fan, which could in principle reduce the size of the triads. 

Nevertheless, we believe that our new method can be used to obtain the vast majority of the 
triads within the given box, by collecting the triads obtained from combinations of small-enough 
primes pj,qj,rj. We provide numerical support of this belief in the next paragraphs. Of course, 
we cannot wait the 15 years in order to check explicitly the percentage of triads in the box 
\k\, \l\ < 5000 that are generated with our new method. For this reason, to prove the point we 
performed a brute-force search to find all resonant triads in the box \k\, \l\ < 100. This took 
about 40 minutes using a desktop computer. In contrast, on the same computer our method 
takes only 1/2 of a second to generate the resonant triads within the box (along with many more 
triads outside the box), using a non-parallel computation of the triads generated by the following 

16 numbers: N = 4 x p 3 -x 7 2 x 5 2 and N = 4 x p 2 x 7 2 x 5 2 , with j = 1, . . . , 8, along with 
the Pythagorean triples generated by u = 2,v = ±1, which is used to generate extra rational 
solutions of the sum of squares problem. 

Following the insight provided by the "small-box" case, we have developed a search method 
that looks at the triads in the box \k\, \l\ < 5000 that are generated by: 

The 1600 numbers iV = 4 x p 3 -x q\ x rf, with j = 1, . . . , 100, k = 1, ... ,4 and / = 
1, . . . ,4, combined appropriately with the 42 Pythagorean triples generated by 1 < \v \ < u < 8. 
Computing time: 287 seconds. Triads obtained: 510 irreducible triads within the box. 
The 627 numbers iV = 4 x p 3 x 7 2 x 5 2 with j = 101, . . . , 727, combined appropriately with 
the 42 Pythagorean triples generated by 1 < \v\ < u < 8. Computing time: 117 seconds. Triads 
obtained: 66 irreducible triads within the box. 

The 1792 numbers N = 4 x p 2 x q\ x rf , with j = 1, . . . , 7, k = 1, . . . , 16 and I = 
1, . . . , 16, combined appropriately with the 42 Pythagorean triples generated by 1 < \v \ < u < 8. 
Computing time: 665 seconds. Triads obtained: 134 irreducible triads within the box. 
The 120 numbers N = 4 x pjp k x 7 2 x 5 2 , with j = 1, 2, 3 and k — j + 1, . . . , j + 40, combined 
appropriately with the 42 Pythagorean triples generated by 1 < \v \ < u < 8. Computing time: 
122 seconds. Triads obtained: 160 irreducible triads within the box. 

We have tested higher-degree combinations of primes and we have obtained up to 20 more 
triads, however for simplicity of presentation we do not consider these here. In total, our search 
takes about 20 minutes on an 8-core desktop computer using Mathematica. We obtain a total 
of 870 irreducible triads, which leads to a total of 6794 triads (including reducible ones) in the 
box \k\, \l\ < 5000. A plot of (i) the number of irreducible resonant triads and (ii) the number 
of total resonant triads, contained in a box as a function of the box size, is shown in Figure [TJ 
left panel. The right panel plots the same variables in log-log scaling. 
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Resonant triads within box of size L Resonant triads within box of size L (Log-Log plot) 




FIG. 1. See text. 



It is apparent from figure [I] right panel, that the total number of resonant triads (including 
reducible ones) as a function of the box size, behaves as a power law of the size, in the form: 
T = CL a , where a « 1.2 and C « 0.263 (this is obtained using fit interval 200 < L < 4000). 
Notice that our numerical method gives an exact estimate for the number of triads for L = 100, 
and the estimation is by construction more accurate for smaller L than for larger L. So, obtaining 
a power law for the number of triads over a range of values of L is a good indication of the 
accuracy of our estimations. 

Another aspect of the set of resonant triads within the box of size L = 5000 is the connectivity 
of the web of resonant triads (including reducible ones). In figure [2] we show a diagrammatic plot 
of the clusters that appear within a box of size L = 200 (for higher values of L the plots become 
prohibitively populated). We notice the predominance of isolated triads. This feature persists for 
larger box sizes, and the isolated triads contribute with approximately 50% of the total number 
of resonant triads. However, there is a nontrivial distribution of cluster sizes within a given box 
size. We show in figure [3] top left panel, the log-log plot of the distribution of number of clusters 
(y-axis) as a function of cluster size (i.e., the number of modes in the cluster), for the collection 
of resonant clusters within the box of size L = 5000. For example, isolated triads have a cluster 
size n mo des = 3 and butterflies have cluster size n modcs = 5. For reference, we plot the size of 
the biggest cluster as a function of box size L in figure [3] top right panel. It is evident that the 
biggest cluster contributes with only a fraction of the modes involved in resonant interactions 
(about 1%). Finally, the total number of modes n tot involved in resonant interactions within a 
given box of size L, as a function of L, behaves in a way that is similar to the total number of 
triads in the box. In fact, due to the fact that resonant triad connections occur via one shared 
mode at a time, it can be shown analytically that the total number of modes involved in resonant 
interactions is bounded between 2n triads + 1 and 3n tria d s , where n triads is the total number of 
triads. The origin of these bounds is as follows: the lower bound comes from the assumption that 
all modes are connected in a single cluster, with connections between triads via a single common 
mode, as the clusters in figure [2] The upper bound is the case when all modes form isolated 
triads. Finally, in figure [3] lower panel, we show the plot of the total number of modes involved 
in resonant interactions n mo d es as a function of box size L, along with the respective bounds. 
Empirically, the relation between n mo des and n tria ds is found to be linear (figure not shown). A 
linear fit (fitting interval 200 < L < 5000) gives a relation n modcs ~ 2.6n triads — 60. 
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FIG. 2. Symbolic plot of the CHM resonant clusters generated within the box of size L = 200. 



V. QUASI-RESONANT TRIADS 

Quasi-resonant triads are defined by the system of equations 

h + k 2 = k 3 (38) 

h + h = h, (39) 

along with the inequality 

\ui +co 2 - co 3 \ < 5, (40) 
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FIG. 3. See text. 



where the bound 5 is known as the detuning level allowed for a triad. 

Quasi-resonant triads play a crucial role in the dynamics of CHM and are more important 
physically than the resonant triads. This is because in any realistic setting the dispersion relation 
has an experimental error, so in fact every physically sensible triad is quasi-resonant. Moreover, 
the amplitude of the wave oscillations is always finite, not infinitesimally small. Therefore the 
time scale of non-resonant triads is not infinitesimally small, and in fact it can be comparable 
with the time scale of the resonant nonlinear oscillations, provided 5 is small enough. The 
accepted conclusion is that a full understanding of the dynamics of a wave system requires the 
understanding of the web of quasi-resonances, rather than just the web of resonances. Notice 
that the web of quasi-resonances for a given dispersion relation is robust (i.e. stable) under small 
perturbations of the dispersion relation, whereas the web of exact resonances is unstable under 
small perturbations. 

These considerations would appear to imply that our method presented in this paper is useless. 
Fortunately, this is incorrect because our method can be used directly to generate quasi-resonant 
triads, in a hierarchical way in the sense that the triads generated have mismatch u\ + u 2 — w 3 
that is small. 



Numerical Method to generate quasi-resonant triads within a given box, starting from 
exact resonant triads of any size. This is based on the fact that the dispersion relation in 
CHM, equation ([3]) is homogeneous under overall re-scaling of the wavenumbers fci, h, k 2 , h, k 3 , l 3 
by any constant (not necessarily integer). Our previous method to generate exact resonant triads 
starts with a number iV and computes all possible representations of the prime expansion of this 
number N. Typically, one gets a lot of triads that are outside the given box. While these triads 
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were discarded in the previous method, the new method uses all resonant triads available. 

From here on we take (3 = — 1 for simplicity in the dispersion relation ([3]). Let (Ki,Li), 
(K 2 , L 2 ), (K s , L 3 ) be an irreducible resonant triad, with Kj,Lj integers. Then the re-scaled 
triad (aKi, aLi), (aK 2 ,aL 2 ), (aK 3 , aL 3 ) is resonant, for any ael. However this triad is not 
necessarily integer, so we need to approximate the scaled wavevectors to nearby integers, keeping 
in mind that equations ([5]) and should be satisfied. This approximation will generate an 
error in the individual frequencies, so equation (|7]) is not satisfied anymore but inequality (40) is 
satisfied, with detuning levels that depend on the accuracy of the approximation from re-scaled 
wavevectors to integer wavevectors. 

In practice, suppose we want to generate quasi-resonant triads within a box of size L. We 
take an irreducible triad that is outside the box and re-scale it to just fit within the box (i.e., the 
maximum modulus of the six re-scaled components is now equal to L). The re-scaled components 
are written as real numbers in, say, double precision. This is our fundamental triad. Next, we 
take re-scaled copies of our fundamental triad, with box norms equal to L — 1, L — 2, . . . , 1. For 
each of these copies we approximate the triad to nearby integer triads, respecting the first two 
resonance conditions and (|6]). The resulting set of triads is formed by quasi-resonant triads 
only, and the higher the box norm the smaller the size of the corresponding detuning level. 

We repeat the process for every irreducible triad outside the box. In this way we can get a 
huge number of quasi-resonant triads (ordered by detuning level) starting from a comparatively 
small set of resonant triads. 

The advantage of this new method is that it provides, in a matter of seconds, a set of quasi- 
resonant triads whose detuning levels are distributed smoothly about the origin. In contrast, the 
brute-force search for quasi-resonant triads takes a long time to generate a smooth distribution, 
of the order of the time it takes to scan all possible triads, resonant or non-resonant (15 years 
for a box of size 10000). 



Numerical Results. We provide a computation of a sample of the quasi-resonant triads within 
the box of size L = 100. Although for this box size the problem can still be treated by brute- 
force search (there are only L 4 = 10 8 triads in total), the choice of a small box size allows us 
to illustrate effectively the main results in terms of connectivity of the quasi-resonant clusters. 
We start generating the irreducible resonant triads generated by the prime expansions of iV = 
4 x Pj x q% x rf, with j — 1, ... ,8, k — 1, . . . , 4 and / = 1, . . . , 4, combined appropriately 
with the 2 Pythagorean triples generated by u — 2, v — ±1. These are 14854 irreducible triads in 
total (computation time = 4 seconds). Some of these triads are within the box of size L = 100, 
but the majority are outside. We pick at random, a small sample of 40 triads out of the 14854 
irreducible triads, so that the sample is symmetric under mirror symmetry k y — > — k y (this is 
in order to allow for eventual connections of clusters with their mirror images). We apply to 
these 40 triads the re-scaling algorithm explained above to generate quasi-resonant triads. As a 
result, 40, 434 quasi-resonant triads within the box are formed. Their distribution in terms of the 
values of frequency mismatch is plotted in figure [4j left panel. The histogram appears uneven 
only because of the convention we use to store the representative triads (0 < k\ < k 2 < k 3 ), 
and physically the distribution must be understood as the symmetrisation of the histogram. The 
corresponding group of clusters obtained is impossible to plot symbolically in a clear way, because 
40434 connected modes are involved. However, we can study quantitatively the change of the 
connectivity properties of the clusters as a function of the allowed detuning level. Consider figure 
[4| right panel (which is analogous to figure [3] lower panel, except that now the independent 
variable is the detuning level 5). For small detuning, below 5 = 0.00004, the connections 
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are still predominantly between isolated triads and connected clusters with one-common-mode 
connectivity. However, the clusters formed when 5 goes beyond about 0.00005, violate the lower 
bound that assumes that connectivity, so we deduce that two-common-mode connectivity is 
starting to prevail. It is easy to show that this new connectivity has a lower bound for the 
number of modes, corresponding to the linear relation n mo d es = ^triads + 2. So, the conclusion is 
that the clusters that form when allowed detuning is high enough, show predominantly the new 
type of triad connectivity. 

Onset of turbulence. A plot, in (k,l) wavevector space, of the modes in the biggest cluster 
obtained when the allowed detuning is 5 = 0.0002, is shown in figure [5] left panel. In the middle 
panel, we plot the biggest cluster when detuning is 5 = 0.00004 and in right panel, we plot 
the biggest cluster when § = 0.00003. It is evident that as the allowed detuning decreases, only 
high wavenumbers in size are allowed to interact and the angular distribution of active modes 
becomes more anisotropic. 

The previous analyses showed that, as the allowed detuning goes beyond some threshold, the 
quasi-resonant triads tend to get connected into one big cluster, with connectivity that shifts 
gradually from a one-common-mode regime to a two-common-mode regime. Let us accept the 
hypothesis that, in a statistically invariant physical system, the typical amplitude of oscillations 
J < ip(x, y, t) 2 > (averaged over x,y,t) activates quasi-resonant triads with allowed detuning 
5 oc y/< ip(x,y,t) 2 >. Evidence in favour of this hypothesis is given by a scaling argument in 
CHM equation ([I]), where the amplitude scales as the inverse of time. Therefore the detuning 
scales as the amplitude. We deduce that if, in a physical system, we increase gradually the 
amplitude of the oscillations (say, via forcing or via manipulating the initial conditions), then 
there is a threshold amplitude below which energy transfers throughout the spectrum of scales 
and directions are not permitted, and above which they are permitted and in fact favoured by 
the fact that the connectivity is via two common modes. 



VI. NEW DEVELOPMENTS: ASPECT RATIO NOT EQUAL TO ONE 

In this section we consider briefly the case when aspect ratio is not equal to 1. In other words 
we will be considering the case where 
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where / = fiyffz, fi & Q and f 2 E N square free. The choice of / so that f 2 is rational is 



important: otherwise the resonance condition analogous to equation (10) would split into two 
or more independent equations, depending on the algebraic degree of /. For example, if / was 
transcendental then all non-zero solutions of the resonance conditions would involve zonal modes 
(e.g., k 3 = 0) which do not interact physically in an exact resonant triad. 

The analysis and transformations for the case f 2 E Q are analogous to the case f — 1, with 
some minor changes. We omit the intermediate steps and just present the final form of the 
mappings. With the assumption k 3 ,ki 0, we define 



D 



k\ k 2 + f 2 l 2 



and we have two cases: 



A. Case D = 



In this case, the ratio must be a pure-cube rational and we obtain a triad of the form 

(h,h), (k2,k), (k 3 ,l 3 ) with 

/ 7 \l/3 

{h, h) = - iji^j (-/%, h), (k 2 , l 2 ) = (k 3 -h u h- h), (41) 
where an overall scaling factor might be needed in order that all wavenumbers be integer. 



B. Case D ^ 

Provided the ratio l 3 / f 2 k 3 is not a pure-cube rational, we can produce a rational bijective 
mapping from an integer triad (k\,li), (^2,^2), (^3,^3) satisfying all resonance conditions, to the 
variables of an elliptic curve 



X 3 -2DX 2 + 2DX-D 2 = fY 2 



(42) 
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The mapping, bijective up to re-scaling of the triad wavenumbers, is given explicitly by: 

x = h kf + f 2 l\ Y = — k ^ 1 - klls D= — kdkl + (43) 
k\ k 2 + f 2 l 2 h\ k 2 + f 2 l\ k\ k 2 + f 2 l 2 



and with inverse 

fci X h X / L> \ Z 3 L> - 1 



k 3 D 2 + f 2 Y 2 ' k 3 f 2 Y V -D 2 + f 2 Y 2 J ' fc 3 / 2 F 
and (k 2 ,h) = (h - ki,l 3 ~ h)- 



(44) 



C. Classification of solutions of triad equations in terms of Fermat's theorem of sums of 

squares (case D / 0) 

We consider the case D^O described in the previous Subsection. 



We can rewrite the elliptic curve (42) as: 

f 2 Y 2 + (D + X 2 -X) 2 = X 2 (X 2 -X + l) , (45) 
and we can divide by X because 1^0, obtaining: 

f (^) 2 +(^ + X - 1 ) 2 = X2 - X + 1 - ( 46 ) 
The last expression is a quadratic form that is best written in diagonal form by defining 

Tfl ~\~ ft 

X = , m,neZ, (47) 

m — n 

so we get the equation: 

f 2 I — +(— - + 2m) =3m 2 + n 2 . 48) 

\ m + n J \ m + n J 

Finally, using the definition of the aspect ratio / = fwffz with /i G Q and / 2 6N square free, 
we obtain the equation: 

\ m+n J \ m + n J 

The problem is thus reduced to finding representations of integers as sums of the form 3m 2 + n 2 
with m,n integers, and of the form f 2 s 2 + q 2 with s,q rationals, where f 2 is a square-free natural 
number. The solution to this problem depends of course on the explicit value of f 2 , but the 
method is straightforward and computable, due to the so-called H a sse- Minkowski theorem, see 
[2"21 pages 61-69]. Essentially, the solution algorithm will be very similar to the one shown in 
Section [Til B[ The details will be presented in a forthcoming paper. 
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Appendix A: Numbers of the form 3m 2 + n 2 



The detailed discussion on this topic is available in [21, Chapter 1] and [221 Section 7]. We 
briefly describe Fermat's approach to finding, for a given integer N, all integers m,n such that 
N = 3m 2 + n 2 . The basic observation is that all we need is to consider such problem when iV 
is a prime number. Then a general N is written as a product of primes and the corresponding 
question can be solved by using the Brahmagupta identity iteratively: 

(Ma 2 + b 2 )(Mc 2 + d 2 ) = M (ad + bc) 2 + (bd -Mac) 2 = M(ad- be) 2 + (bd + M ac) 2 , (Al) 

valid for any square-free M e Z and a, b,c,de Q. 

Fermat's argument shows that a prime P can be written in the form P = 3m 2 + n 2 , if and 
only if P — 3 or P — 3k + 1, where k is some integer. In such case the coprime integers m,n 
are unique up to signs. Similarly, any positive integer power of such primes has a unique (up to 
sign) representation in the form 3 m 2 + n 2 . For example, 7 = 3 x l 2 + 2 2 , 7 2 = 3 x 4 2 + l 2 and 
7 3 = 3 x 9 2 + 10 2 in unique ways. All these can be obtained from application of the Brahmagupta 
identity iteratively starting from the expansion for 7 = 3 x l 2 + 2 2 . Non-coprime integers are 
obtained trivially from the coprime solutions. For example, we can also write 7 3 = 7 2 x 7 = 
7 2 x (3 x l 2 + 2 2 ) = 3 x 7 2 + 14 2 . 

By allowing for the special case 4 = 3 x l 2 + l 2 , one can state the result below due to Fermat 
[2U Chapter 1]: 

Theorem 1. A general positive integer N is expressible in the form A" = 3 m 2 + n 2 , if and only 
if the prime expansion of A" contains no odd powers of primes of the form P = 3 k — 1 , with k 
integer. It can contain odd powers of 4, 3 and primes P — 3 k + 1. Moreover, the ways in which 
the integer A" can be represented as N = 3m 2 + n 2 with m,n integers, can be computed by 



repeated iteration of the Brahmagupta identity (Al) to each product of primes. 



For example, the number A^ = 4 x 3 x 7 = 84 can be written as 84 = 3x5 2 + 3 2 = 3x l 2 + 9 2 
But neither AT = 5 nor A^ = 11 can be written in such form. 



Remark. For our computational purposes, for a given positive integer A" according to Theorem 
1, it is important to have an algorithm to compute all possible values of the integers m',n' in 
A" = 3 [m') 2 + (n') 2 . Such algorithm is based on the complex representation of prime factors 
P = 3 m 2 + n 2 = (n + \/3im)(n — \f3im) , where i = and is equivalent to the iterative 

application of Brahmagupta identity QAip . 



Appendix B: Numbers of the form s 2 + q 2 



This is the best known problem of this type, leading to the so-called Pythagorean triples, see 
[23 chapter V] and [20, Chapter VI]. The relevant questions are: 
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Which positive integers N can be represented in the form N = s 2 + q 2 , with s, q rationals? In 
how many ways can we write such representation? 

In this case s,q can be rationals because there are no restrictions on Y and D in equation 



(24). However, the question in terms of rationals can be reduced to a question in terms of 



integers, by writing the number 1 as a sum of squares of two rationals: 

Wu, v G Z. 



> \ - / 2 2\ 2 

2uv \ ( u • — v ' x 



U 2 + V 2 J \U 2 + V 2 

This equation is easy to verify (directly or via Brahmagupta identity) and it leads to the 
Pythagorean triples (2 uv,u 2 — v 2 , u 2 + v 2 ), which are the integer sides of a right triangle. Us- 



ing equation (29) it is possible to write down a partial solution to the question above, in the form: 



Theorem 2. Let N = s 2 + q 2 with N G Z integer and s, q G Q. Then, there exist u, v G Z and 
S, Q G Z such that 

2uvQ + (u 2 -v 2 )S -2uvS + (u 2 - v 2 ) Q 

S= 2~^ 2 ' V = 2~s 2 ' B1 

u z _|_ v z u z _|_ y z 

and N = S 2 + Q 2 . 

Now that the problem has been reduced to finding the integer representations of iV in the 
form N = S 2 + Q 2 , we can use Fermat's results: 

Theorem 3. [23, Section I, Chapter V] A general positive integer iV is expressible in the form 
N = S 2 + Q 2 , if and only if the prime expansion of iV contains no odd powers of primes of the 
form P = Ak — 1 , with k integer. It can contain odd powers of 2 and primes P = Ak + 1. 
Moreover, the ways in which the integer N can be represented as N = S 2 + Q 2 with S, Q 



integers, can be computed by repeated iteration of the Brahmagupta identity (Al) to each 
product of primes. 

Numbers of the form 3m 2 + n 2 and s 2 + q 2 . Combining the results encapsulated in Theorems 
1, 2 and 3, we obtain the following 

Corollary 4. A general positive integer iV is expressible in the forms iV = 3 m 2 + n 2 = s 2 + q 2 
with m,n G Z and s, q G Q, if and only if the only odd powers of primes appearing in the 
prime expansion of N are those of the primes of the form P = 12k + 1, where k is some 
integer. As a special case, odd powers of 4 are also allowed in the prime expansion of N. 
Moreover, the different solutions for m, n and s,q can be fully classified in terms of the basic 



representations of the prime factors in N and the Pythagorean triples appearing in equation (29) 



-2uvS+(u 2 -v 2 )Q ^ 2 
u 2 _^_ v 2 



For example, the integer N = 13 is expressible as iV = 3 x 2 2 + l 2 = ( 2uvQ $£? P)S ) + 



, with (S,Q) = (3,2) and u, v arbitrary integers. 
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